The field of the invention is medical imaging and particularly, methods for reconstructing images from acquired image data.
Magnetic resonance imaging (“MRI”) uses the nuclear magnetic resonance (“NMR”) phenomenon to produce images. When a substance such as human tissue is subjected to a uniform magnetic field (polarizing field B0), the individual magnetic moments of the nuclei in the tissue attempt to align with this polarizing field, but precess about it in random order at their characteristic Larmor frequency. If the substance, or tissue, is subjected to a magnetic field (excitation field B1) that is in the x-y plane and that is near the Larmor frequency, the net aligned moment, Mz, may be rotated, or “tipped,” into the x-y plane to produce a net transverse magnetic moment Mxy. A signal is emitted by the excited nuclei or “spins,” after the excitation signal B1 is terminated, and this signal may be received and processed to form an image.
When utilizing these “MR” signals to produce images, magnetic field gradients (Gx, Gy, and Gz) are employed. Typically, the region to be imaged is scanned by a sequence of measurement cycles in which these gradients vary according to the particular localization method being used. The resulting set of received MR signals are digitized and processed to reconstruct the image using one of many well known reconstruction techniques.
The measurement cycle used to acquire each MR signal is performed under the direction of a pulse sequence produced by a pulse sequencer. Clinically available MRI systems store a library of such pulse sequences that can be prescribed to meet the needs of many different clinical applications. Research MRI systems include a library of clinically-proven pulse sequences and they also enable the development of new pulse sequences.
The MR signals acquired with an MRI system are signal samples of the subject of the examination in Fourier space, or what is often referred to in the art as “k-space.” Each MR measurement cycle, or pulse sequence, typically samples a portion of k-space along a sampling trajectory characteristic of that pulse sequence. Most pulse sequences sample k-space in a raster scan-like pattern sometimes referred to as a “spin-warp,” a “Fourier,” a “rectilinear,” or a “Cartesian” scan. The spin-warp scan technique employs a variable amplitude phase encoding magnetic field gradient pulse prior to the acquisition of MR spin-echo signals to phase encode spatial information in the direction of this gradient. In a two-dimensional implementation (“2DFT”), for example, spatial information is encoded in one direction by applying a phase encoding gradient, Gy, along that direction, and then a spin-echo signal is acquired in the presence of a readout magnetic field gradient, Gx, in a direction orthogonal to the phase encoding direction. The readout gradient present during the spin-echo acquisition encodes spatial information in the orthogonal direction. In a typical 2DFT pulse sequence, the magnitude of the phase encoding gradient pulse, Gy, is incremented, ΔGy, in the sequence of measurement cycles, or “views” that are acquired during the scan to produce a set of k-space MR data from which an entire image can be reconstructed.
There are many other k-space sampling patterns used by MRI systems. These include “radial”, or “projection reconstruction” scans in which k-space is sampled as a set of radial sampling trajectories extending from the center of k-space. The pulse sequences for a radial scan are characterized by the lack of a phase encoding gradient and the presence of a readout gradient that changes direction from one pulse sequence view to the next. There are also many k-space sampling methods that are closely related to the radial scan and that sample along a curved k-space sampling trajectory rather than the straight line radial trajectory.
An image is reconstructed from the acquired k-space data by transforming the k-space data set to an image space data set. There are many different methods for performing this task and the method used is often determined by the technique used to acquire the k-space data. With a Cartesian grid of k-space data that results from a 2D or 3D spin-warp acquisition, for example, the most common reconstruction method used is an inverse Fourier transformation (“2DFT” or “3DFT”) along each of the 2 or 3 axes of the data set. With a radial k-space data set and its variations, the most common reconstruction method includes “regridding” the k-space samples to create a Cartesian grid of k-space samples and then performing a 2DFT or 3DFT on the regridded k-space data set. In the alternative, a radial k-space data set can also be transformed to Radon space by performing a 1DFT of each radial projection view and then transforming the Radon space data set to image space by performing a filtered backprojection.
Depending on the technique used, many MR scans currently used to produce medical images require many minutes to acquire the necessary data. The reduction of this scan time is an important consideration, since reduced scan time increases patient throughout, improves patient comfort, and improves image quality by reducing motion artifacts. Many different strategies have been developed to shorten the scan time.
In an x-ray computed tomography (“CT”) system, an x-ray source projects a fan-shaped beam that is collimated to lie within an x-y plane of a Cartesian coordinate system, termed the “image plane.” The x-ray beam passes through the object being imaged, such as a medical patient, and impinges upon an array of radiation detectors. The intensity of the transmitted radiation is dependent upon the attenuation of the x-ray beam by the object and each detector produces a separate electrical signal that is a measurement of the beam attenuation. The attenuation measurements from all the detectors are acquired separately to produce what is called the “transmission profile,” or “attenuation profile,” or “projection.”
The source and detector array in a conventional CT system are rotated on a gantry within the imaging plane and around the object so that the angle at which the x-ray beam intersects the object constantly changes. The transmission profile from the detector array at a given angle is referred to as a “view” and a “scan” of the object comprises a set of views made at different angular orientations during one revolution of the x-ray source and detector. In a 2D scan, data is processed to construct an image that corresponds to a two dimensional slice taken through the object. The prevailing method for reconstructing an image from 2D data is referred to in the art as the filtered backprojection technique. This image reconstruction process converts the attenuation measurements acquired during a scan into integers called “CT numbers,” or “Hounsfield units,” which are used to control the brightness of a corresponding pixel on a display.
According to standard image reconstruction theories, in order to reconstruct an image without aliasing artifacts, the sampling rate employed to acquire image data must satisfy the so-called Nyquist criterion, which is set forth in the Nyquist-Shannon sampling theorem. Additionally, in standard image reconstruction theories, no specific prior information about the image is needed. On the other hand, when some prior information about the desired, or target, image is available and appropriately incorporated into the image reconstruction procedure, an image can be accurately reconstructed even if the Nyquist criterion is violated. For example, in parallel-beam CT imaging, if one knows that a desired target image is circularly symmetric and spatially uniform, only one projection view of parallel-beam projections is necessary to accurately reconstruct the linear attenuation coefficient of the object. As another example, if one knows that a desired target image consists of only a single point, then only two orthogonal projections that intersect at that point are needed to accurately reconstruct the image point.
Thus, if prior information is known about the desired target image, such as if the desired target image is a set of sparsely distributed points, it can be reconstructed from a set of data that was acquired in a manner that does not satisfy the Nyquist criterion. Put more generally, knowledge about the sparsity of the desired target image can be employed to relax the Nyquist criterion; however, it is a highly nontrivial task to generalize these arguments to formulate a rigorous image reconstruction theory.
The Nyquist criterion serves as one of the paramount foundations of the field of information science. However, it also plays a pivotal role in modern medical imaging modalities such as MRI and x-ray CT. When the number of data samples acquired by an imaging system is less than the requirement imposed by the Nyquist criterion, artifacts appear in the reconstructed images. In general, such image artifacts include aliasing and streaking artifacts. In practice, the Nyquist criterion is often violated, whether intentionally or through unavoidable circumstances. For example, in order to shorten the data acquisition time in a time-resolved MR angiography study, undersampled projection reconstruction, or radial, acquisition methods are often intentionally introduced.
The common problem presented in imaging is of estimating a discrete approximation of an image that gave rise to the observed measurements. If Φ denotes a linear sensing operator, such as the partial discrete Fourier Transform (“DFT”) matrix common in magnetic resonance imaging, this estimation problem becomes one of recovering a discrete signal, s, from a corrupted measurement vector:f=Φs+n  Eqn. (1);
where n represents the statistics of a noise process and f is the measurement vector. In MRI, the measurement vector is representative of acquired k-space data, whereas in CT, the measurement vector is representative of the acquired x-ray projection data. It is noted that generic a priori knowledge about the signal class, of which the discrete signal, s, is a member, may be considered. It can be assumed that n is drawn from an independent and identically distributed (“i.i.d.”) Gaussian processes.
Given Eqn. (1), the maximum likelihood (“ML”) estimator of the discrete signal, s, is given by the solution to the ordinary least squares (“OLS”) problem:
                                          u            ML                    =                                    arg              ⁢                                                          ⁢                                                min                  u                                ⁢                                                                                                                        Φ                        ⁢                                                                                                  ⁢                        u                                            -                      f                                                                            2                  2                                                      =                                                            (                                      Φ                    *                    Φ                                    )                                                  -                  1                                            ⁢              Φ              *              f                                      ;                            Eqn        .                                  ⁢                  (          2          )                    
where u is a desired estimate of the discrete signal, s. Note that the uniqueness of uML is only guaranteed if Φ is full rank. In fact, even if Φ is full rank, it may very well be ill-conditioned. To counter potential ill-posedness, the OLS problem is often stabilized by imposing constraints based on a priori knowledge or assumptions about the discrete signal, s.
Sparsity is a very “powerful” form of a priori knowledge. An image, I, is said to be sparse if, given some “sparsifying” transform, Ψ, such as a wavelet or spatial gradient transform, S|Ω|, where S=supp(ΨI) and Ω=dom(I). Given a limited set of potentially noisy measurements, the search for the sparsest image possessing a prescribed degree of fidelity to the measurement set is defined by the following l0-minimization problem:
                                                        min              u                        ⁢                                                                                                  Ψ                    ⁢                                                                                  ⁢                    u                                                                    0                            ⁢                                                          ⁢                              s                .                t                .                                                                  ⁢                                                                                                                        Φ                        ⁢                                                                                                  ⁢                        u                                            -                      f                                                                            2                  2                                                              ≤                      ɛ            2                          ;                            Eqn        .                                  ⁢                  (          3          )                    
where ε2 is the variance of the noise, n, and ∥v∥0 indicates the l0-norm of a vector, v, which has the form:
                                                      v                                0                =                              ∑                          x              ∈              Ω                                ⁢                      1            ⁢                                          (                                                                                                v                      ⁡                                              (                        x                        )                                                                                                  >                  0                                )                            .                                                          Eqn        .                                  ⁢                  (          4          )                    
Thus, the l0-norm of a vector, v, provides a count of the non-zero terms in the vector, v.
It has recently been shown that, for the specific case in Eqn. (3) where the sparsifying transform, Ψ, is set to the identity matrix, and ε=0, only O(S log(|Ω|)) random samples are necessary to guarantee with significant probability that the solution of Eqn. (3) is exactly s. While these sampling conditions are theoretically alluring, Eqn. (3) requisites combinatorial optimization (Eqn. (3) is, technically speaking, NP-complete) and, thus, is of little practical use in modern imaging.
Recently, a new mathematical framework for image reconstruction termed “compressed sensing” (“CS”) was formulated. In compressed sensing, only a small set of linear projections of a sparse image are required to reconstruct a quality image. The theory of CS is described by E. Candès, J. Romberg, and T. Tao, in “Robust uncertainty principles: Exact Signal Reconstruction from Highly Incomplete Frequency Information,” IEEE Transactions on Information Theory 2006; 52:489-509, and by D. Donoho in “Compressed Sensing,” IEEE Transactions on Information Theory 2006; 52:1289-1306, and is disclosed, for example, in U.S. Pat. No. 7,646,924.
Consider the following generalization of Eqn. (3):
                                                        min              u                        ⁢                                          P                ⁡                                  (                                      Ψ                    ⁢                                                                                  ⁢                    u                                    )                                            ⁢                                                          ⁢                              s                .                t                .                                                                  ⁢                                                                                                                        Φ                        ⁢                                                                                                  ⁢                        u                                            -                      f                                                                            2                  2                                                              ≤                      ɛ            2                          ;                            Eqn        .                                  ⁢                  (          5          )                    
where P( . . . ) is a separable metric prior functional. The work of Candès, et al., and Donoho produced both exact and approximate recovery sampling conditions for P(v)=∥v∥1, a convex functional, similar to those associated with Eqn. (3), albeit with a higher sampling constant. This area of study is now generally known as compressed, or compressive, sensing.
In aspiration of reducing sampling requirements closer to the l0-associated theoretical limit, quasi-convex, almost everywhere (“a.e.”) differentiable prior functionals, such as P(v)=∥v∥pp with 0<p≦1 have been investigated recently. Such prior functionals more closely resemble the target l0-norm than does l1; however, while these methods have the potential to outperform analogous convex techniques, in practice, they cannot guarantee numerical achievement of global minima since they are non-convex in nature.
The field of medical imaging can benefit from compressed sensing. In x-ray computed tomography (“CT”), fewer measurements translates to a lower radiation dose received by the patient. Similarly, in magnetic resonance imaging (“MRI”), decreasing the required number of measurements necessary to form an image allows for faster exams, which improves patient comfort and thereby minimizes the likelihood of subject motion. Additionally, clinical throughput is increased which potentially translates to reduced patient costs.
Since its inception, substantial effort has been made to decrease the required scan time in MRI. Early techniques, such as echo-planar imaging (“EPI”), traverse the entirety of k-space during a single repetition time (“TR”) period, thereby offering a dramatic increase in speed, albeit at the expense of demanding hardware performance and significantly lowered signal-to-noise (“SNR”) levels.
Alternatively, if the signal of interest is assumed to be strictly real, Hermitian symmetry of k-space can be exploited such that it is only necessary to perform measurements over half of the spectral grid in k-space. In practice, however, the strict reality assumption is violated due to magnetic field inhomogeneities; chemical shift and magnetic susceptibility effects; as well as artifacts from physiological motion and flow dynamics. As a result of these effects, phase errors must be corrected prior to making the Hermitian assumption. If image phase is assumed to be smoothly varying over space, the standard approach to phase correction involves measuring a symmetric low-frequency band of k-space and estimating the image space phase solely from this restricted measurement set. The solution magnitude is then derived from a moderately undersampled subset of k-space which also includes the small support used in the phase estimation step. A complex image is then formed by conjoining the image magnitude and phase estimates, with only the real portion of this generated image being retained as the solution. Due to the necessity of acquiring a symmetric, low-frequency spectral band for use in the phase estimation process, methods which rely on Hermitian symmetry such as projections onto convex sets (“POCS”) and homodyne detection can only decrease the number of required measurements by less than half of that delimited by Shannon's theorem.
In both MR and CT imaging, any reduction in scan time offers potential benefits including high temporal resolution for observing physiologic processes, reduction of motion artifacts, lower radiation dose (for CT), and improved patient comfort and throughput. With the advent of compressive sensing techniques, certain classes of images can be accurately reconstructed from highly undersampled, and thus quickly acquired, data by solving a convex l1-minimization problem. For example, following the recent development of compressed sensing, it was shown that clinical MR images possessing a sparse representation in some transform domain can be accurately reconstructed, even when sampled at rates well below the Nyquist limit. This accurate reconstruction was achieved by casting the image reconstruction as a convex l1-minimization problem. While l1-based techniques offer an advantage over Nyquist-limited methods, they nonetheless require a significant degree of oversampling above the theoretical minimum sampling rate allowed by compressed sensing in order to guarantee the achievability of exact reconstruction. In addition, such techniques are computationally intensive.
It would therefore be desirable to provide a method for image reconstruction, in which a greater degree of undersampling than achievable with current l1-based reconstruction techniques is employed without resulting in a decrease in overall image quality.